{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"Script for plotting various distributions of Milky Way clusters.\n",
    "Data sources:\n",
    "Improving the open cluster census. II. An all-sky cluster catalogue with Gaia DR3\n",
    "https://ui.adsabs.harvard.edu/abs/2023A%26A...673A.114H/abstract\n",
    "VizieR Online Data Catalog: Improving the open cluster census. II. (Hunt+, 2023)\n",
    "Data source: https://cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/673/A114\n",
    "Last modification: 23-May-2023\n",
    "See https://cdsarc.cds.unistra.fr/ftp/J/A+A/673/A114/ReadMe for byte-by-byte\n",
    "description of clusters.dat\n",
    "\"\"\"\n",
    "\n",
    "\n",
    "import os\n",
    "\n",
    "from scour import scour\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.ticker import AutoMinorLocator, ScalarFormatter\n",
    "import pandas as pd\n",
    "\n",
    "\n",
    "def optimize_svg(tmp_path, path):\n",
    "    \"\"\"Optimize svg file using scour\"\"\"\n",
    "    with open(tmp_path, \"rb\") as inputfile, open(path, \"wb\") as outputfile:\n",
    "        options = scour.generateDefaultOptions()\n",
    "        options.enable_viewboxing = True\n",
    "        options.strip_comments = True\n",
    "        options.strip_ids = True\n",
    "        options.remove_metadata = True\n",
    "        options.shorten_ids = True\n",
    "        options.indent_type = \"none\"\n",
    "        options.newlines = False\n",
    "        scour.start(options, inputfile, outputfile)\n",
    "\n",
    "\n",
    "PLOTS_DIR = \"../../../plots/stars/\"\n",
    "DATA_DIR = \"../../../data/hunt2023/\"\n",
    "MS = 3\n",
    "UP = 0.1\n",
    "DATA = np.genfromtxt(\n",
    "    DATA_DIR + \"clusters.dat\",\n",
    "    delimiter=(\n",
    "        21, 5, 247, 2, 12,\n",
    "        7, 12, 6, 13, 13,\n",
    "        13, 12, 12, 12, 12,\n",
    "        12, 14, 14, 14, 14,\n",
    "        14, 12, 11, 13, 12,\n",
    "        11, 13, 12, 11, 16,\n",
    "        16, 17, 6, 2, 17,\n",
    "        17, 17, 14, 14, 14, 5,\n",
    "        10, 10, 11, 11,\n",
    "        11, 4, 11, 11,\n",
    "        12, 10, 11, 11, 10,\n",
    "        11, 11, 12, 12, 12,\n",
    "        3, 2, 2, 3, 4),\n",
    "    dtype=[(\"Name\", \"U20\"), (\"ID\", \"int\"), (\"AllNames\", \"U246\"), (\"Type\", \"U1\"), (\"CST\", \"f8\"),\n",
    "        (\"N\", \"int\"), (\"CSTt\", \"f8\"), (\"Nt\", \"int\"), (\"RAdeg\", \"f8\"), (\"DEdeg\", \"f8\"),\n",
    "        (\"GLON\", \"f8\"), (\"GLAT\", \"f8\"), (\"r50\", \"f8\"), (\"rc\", \"f8\"), (\"rt\", \"f8\"), # GLAT E11.4\n",
    "        (\"rtot\", \"f8\"), (\"r50pc\", \"f8\"), (\"rcpc\", \"f8\"), (\"rtpc\", \"f8\"), (\"rtotpc\", \"f8\"),\n",
    "        (\"pmRA\", \"f8\"), (\"s_pmRA\", \"f8\"), (\"e_pmRA\", \"f8\"), (\"pmDE\", \"f8\"), (\"s_pmDE\", \"f8\"),\n",
    "        (\"e_pmDE\", \"f8\"), (\"Plx\", \"f8\"), (\"s_Plx\", \"f8\"), (\"e_Plx\", \"f8\"), (\"dist16\", \"f8\"),\n",
    "        (\"dist50\", \"f8\"), (\"dist84\", \"f8\"), (\"Ndist\", \"int\"), (\"globalPlx\", \"int\"), (\"X\", \"f8\"),\n",
    "        (\"Y\", \"f8\"), (\"Z\", \"f8\"), (\"RV\", \"f8\"), (\"s_RV\", \"f8\"), (\"e_RV\", \"f8\"), (\"n_RV\", \"int\"),\n",
    "        (\"CMDCl2.5\", \"f8\"), (\"CMDCl16\", \"f8\"), (\"CMDCl50\", \"f8\"), (\"CMDCl84\", \"f8\"),\n",
    "        (\"CMDCl97.5\", \"f8\"), (\"CMDClHuman\", \"U3\"), (\"logAge16\", \"f8\"), (\"logAge50\", \"f8\"),\n",
    "        (\"logAge84\", \"f8\"), (\"AV16\", \"f8\"), (\"AV50\", \"f8\"), (\"AV84\", \"f8\"), (\"diffAV16\", \"f8\"),\n",
    "        (\"diffAV50\", \"f8\"), (\"diffAV84\", \"f8\"), (\"MOD16\", \"f8\"), (\"MOD50\", \"f8\"), (\"MOD84\", \"f8\"),\n",
    "        (\"minClSize\", \"int\"), (\"isMerged\", \"int\"), (\"isGMMMemb\", \"int\"), (\"NXmatches\", \"int\"),\n",
    "        (\"XmatchType\", \"U3\")\n",
    "        ])\n",
    "\n",
    "df = pd.DataFrame(DATA)\n",
    "df[\"Name\"] = df[\"Name\"].str.strip()\n",
    "df[\"age\"] = 10**df[\"logAge50\"] / 10**9\n",
    "df_o = df[df[\"Type\"] == \"o\"]\n",
    "df_m = df[df[\"Type\"] == \"m\"]\n",
    "df_g = df[df[\"Type\"] == \"g\"]\n",
    "ple = df[df[\"Name\"] == \"Melotte_22\"]\n",
    "hya = df[df[\"Name\"] == \"Melotte_25\"]\n",
    "praesepe = df[df[\"Name\"] == \"NGC_2632\"]\n",
    "ruprecht147 = df[df[\"Name\"] == \"Ruprecht_147\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Open clusters:\n",
      "          Name       dist50       age  logAge50   logAge84\n",
      "6439   UBC_584  1680.151735  0.002900  6.462370   6.500426\n",
      "2472   HSC_816   873.315920  0.002937  6.467832   6.545616\n",
      "...        ...          ...       ...       ...        ...\n",
      "3069  HSC_1586  2182.155806  9.030243  9.955699  10.000000\n",
      "1874   HSC_101  2500.904461  9.516442  9.978475  10.000000\n",
      "\n",
      "[6818 rows x 5 columns]\n",
      "Moving groups:\n",
      "          Name        dist50       age  logAge50  logAge84\n",
      "3234  HSC_1795  35988.420830  0.004028  6.605128  7.123780\n",
      "3266  HSC_1837  37103.342194  0.004727  6.674572  7.236075\n",
      "...        ...           ...       ...       ...       ...\n",
      "2249   HSC_565   3057.759936  6.595804  9.819268  9.999873\n",
      "3296  HSC_1868   3202.339059  6.890211  9.838233  9.999987\n",
      "\n",
      "[228 rows x 5 columns]\n",
      "Globular clusters:\n",
      "          Name        dist50       age  logAge50  logAge84\n",
      "4633  NGC_4147  41818.184150  0.004049  6.607337  6.969244\n",
      "4659  NGC_5466  17673.274513  0.010659  7.027730  7.418416\n",
      "...        ...           ...       ...       ...       ...\n",
      "4819  NGC_6838   3983.209664  7.114310  9.852133  9.976394\n",
      "4686  NGC_6121   1771.116638  7.684014  9.885588  9.982500\n",
      "\n",
      "[121 rows x 5 columns]\n"
     ]
    }
   ],
   "source": [
    "with pd.option_context(\"display.max_rows\", 4):\n",
    "    print(\"Open clusters:\")\n",
    "    print(df_o[[\"Name\", \"dist50\", \"age\", \"logAge50\", \"logAge84\"]].sort_values(by=\"logAge50\"))\n",
    "    print(\"Moving groups:\")\n",
    "    print(df_m[[\"Name\", \"dist50\", \"age\", \"logAge50\", \"logAge84\"]].sort_values(by=\"logAge50\"))\n",
    "    print(\"Globular clusters:\")\n",
    "    print(df_g[[\"Name\", \"dist50\", \"age\", \"logAge50\", \"logAge84\"]].sort_values(by=\"logAge50\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Nearest clusters:\n",
      "                Name     dist50       age  logAge50\n",
      "2549         HSC_906  25.084607  0.277234  8.442846\n",
      "1614        FSR_1017  31.180782  0.015689  7.195601\n",
      "4424      Melotte_25  47.190660  0.576836  8.761053\n",
      "4074        HSC_2846  49.965042  0.053457  7.728006\n",
      "7158  beta_Tuc_Group  52.846137  0.031398  7.496899\n",
      "3319        HSC_1900  57.786463  0.026673  7.426064\n",
      "2746        HSC_1152  57.986448  0.130782  8.116549\n",
      "2151         HSC_453  64.796420  0.244329  8.387975\n",
      "3054        HSC_1566  72.347783  0.115044  8.060864\n",
      "2978        HSC_1470  72.655661  0.447163  8.650465\n",
      "3722        HSC_2399  77.190273  0.159707  8.203323\n"
     ]
    }
   ],
   "source": [
    "df_nearest = df[df[\"dist50\"] < 80].sort_values(by=\"dist50\")\n",
    "print(\"Nearest clusters:\")\n",
    "print(df_nearest[[\"Name\", \"dist50\", \"age\", \"logAge50\"]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       ""
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(16, 9))\n",
    "plt.subplots_adjust(left=0.051, bottom=0.06, right=0.98, top=0.955)\n",
    "sct = plt.scatter(df.dist50, df.age, s=MS, label=\"Open clusters\")\n",
    "sct_m = plt.scatter(df_m.dist50, df_m.age, s=MS + 1, c=\"r\", label=\"Moving groups\")\n",
    "sct_g = plt.scatter(df_g.dist50, df_g.age, s=MS + 1, c=\"k\", label=\"Globular clusters\")\n",
    "plt.plot(ple.dist50, ple.age, \"ok\", ms=5)\n",
    "plt.plot(hya.dist50, hya.age, \"og\", ms=5)\n",
    "plt.plot(praesepe.dist50, praesepe.age, \"o\", ms=5, c=\"midnightblue\")\n",
    "plt.plot(ruprecht147.dist50, ruprecht147.age, \"o\", c=\"orange\", ms=5)\n",
    "plt.plot((10, 10000), (4.6, 4.6), \"--y\", lw=2, label=\"The age of the Sun\")\n",
    "\n",
    "plt.text(hya.dist50.iloc[0], hya.age.iloc[0] + UP, \"Hyades\", fontsize=10)\n",
    "plt.text(ple.dist50.iloc[0], ple.age.iloc[0] + UP, \"Pleiades\", fontsize=10)\n",
    "plt.text(praesepe.dist50.iloc[0], praesepe.age.iloc[0] + UP, \"Praesepe\", fontsize=10)\n",
    "plt.text(ruprecht147.dist50.iloc[0] + 2, ruprecht147.age.iloc[0] + UP, \"Ruprecht 147\", fontsize=10)\n",
    "\n",
    "plt.xscale(\"log\")\n",
    "plt.xlim(10, 10000)\n",
    "plt.ylim(-0.02, 9.518)\n",
    "ax.yaxis.set_minor_locator(AutoMinorLocator())\n",
    "plt.title(\"All known clusters in distance-age space (Hunt+, 2023)\")\n",
    "plt.xlabel(\"Distance (pc)\", labelpad=3)  # , fontsize=12\n",
    "plt.ylabel(\"Age (Gyr)\", labelpad=8)\n",
    "plt.legend(loc=\"upper left\")\n",
    "plt.grid(axis=\"both\", which=\"major\", linestyle=\":\")\n",
    "ax.xaxis.set_major_formatter(ScalarFormatter())\n",
    "\n",
    "FILE_EXT = \"png\"\n",
    "PLT_PTH = PLOTS_DIR + \"clusters-dist-age-omg-annotated\"\n",
    "tmp_pth = f\"{PLT_PTH}_.{FILE_EXT}\"\n",
    "pth = f\"{PLT_PTH}.{FILE_EXT}\"\n",
    "plt.savefig(tmp_pth, dpi=120)\n",
    "\n",
    "if FILE_EXT == \"svg\":\n",
    "    optimize_svg(tmp_pth, pth)\n",
    "    os.remove(tmp_pth)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.3"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
